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Abstract 

In the paper we present a functional-discrete method for solving Sturm- 
Liouville problems with potential including function from Li(0, 1) and <5-function. 
For both, linear and nonlinear cases the sufficient conditions providing superex- 
ponential convergence rate of the method are obtained. The question of possible 
software implementation of the method is discussed in detail. The theoretical 
results are successfully confirmed by the numerical example included in the pa- 
per. 
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1 Introduction 

The functional-discrete method (FD-method) was first proposed for the Sturm-Liouville 
problem in [13] . The idea of this approach consists of embedding an original problem 
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into a parametric set of problems with respect to a parameter r in such a way that 
for r = we have a linear eigenvalue problem with piecewise constant coefficients 
and for r = 1 we have the given eigenvalue problem. The transition from t = to 
t = 1 using Taylor series results in a recursive algorithm. Thus, we represent the 
exact solution to the given problem, a pair (eigenvalue, eigenfunction) , as two series, 
the first one — for the eigenvalue, the other one — for the eigenfunction. Then, an 
approximate solution is a pair of the corresponding truncated series. We find the first 
terms for both series via the coefficient approximation method (CAM) using piece- wise 
constant approximations. The CAM was first substantiated by Kryloff and Bogoli- 



aubov ("tronson's" method) in 12 and then developed in [8j[T0j[23j|24] . As a result 
of piece-wise constant approximation, we obtain an "unperturbed" problem which is 
also referenced to as the basic problem. The eigensolution to this problem gives us 
the first terms for the series representations mentioned above. Then, each successive 
term in the series for eigenfunction can be found as the exact solution to a linear 
boundary-value (not eigenvalue!) problem with a scalar parameter. A single value of 
this parameter which provides the solvability of the corresponding BVP gives us the 
next term in the series representation for eigenvalue. This approach allows us to find 
a numerical solution to the given eigenvalue problem with any desired accuracy, that 
is, we can improve the accuracy just by carrying out a few more iterations. 

This technique was developed for the Sturm-Liouville problems with continuous 
potential in [4 
potential in 
space L\ in 



19 



16 



14 



21) 



26 



18 p21 , for the transmission eigenvalue problems with continuous 



25], for the linear eigenvalue problems with potential belonging to 



for the nonlinear eigenvalue problems with continuous potential 
7 p5p7] , For all these cases the sufficient conditions providing the exponential con- 
vergence rate of the method have been found. It was also shown that the convergence 
improves along with the increase of the index of a trial eigenvalue. 

In this article we extend the FD-approach to the case of a discontinuous potential 
which consists of the Delta- function and a function from space L\. Such problems 
are of great interest in [TJ[2[[27}[2S] . 

The paper is organized as follows. In section [2] we state the rigorous formulation 
of the problem we deal with and briefly discuss several known properties of Sturm- 
Liouville problems with distribution potentials. In section|3]the general description of 
the FD-method's algorithm is presented. The main theoretical results are stated and 
proved in section [4j Section [5] is devoted to the question of the software implementa- 
tion for the proposed method. Numerical example is given in section [6j followed by 
section [7] containing the conclusions. 

2 Formulation of the problem 

In the paper we consider the following Sturm-Liouville problerrj^] 

d 2 

— u(x)-{P5(x-a) + q(x)]u(x) + \u(x)-N(u(x))=0, x€(0,l)\{a}, (2.1) 
dx z 

u(0) = u(l) = 0; u'(0) = 1, 



1 In this problem condition u'(0) = 1 can be substituted by f (11(1)) dx = const and the FD- 



method's algorithm presented below can be modified for treating problems of such type. 
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where a G (0, 1), (5 > 0, q(x) G Li [0, 1] , iV(«) = £ a p u p , Vu G R. Here <5(x) denotes 

p=i 

Dirac <5-function or impulse symbol, which can be viewed as the derivative of the 
Heaviside step function (see [H]) 



dx 



H(x)=5(x), H{x) 




(2.2) 



Using (2.2 1 it is easy to verify that the solution u{x) to problem (2.1 1 satisfies the 
following equalities 



u'(x) = [l( x ) u ( x ) ~ Att (as) + N(u(x))] dx, a; G [0, a), 

o 

" 



(2.3) 



Equalities (2.3) imply that function u'(x) is discontinuous at the point x = a: 
u'(a + 0) - u'(a - 0) = /3u(a); 



(2.4) 



Using formulas (2.3) and (2.4) we can conclude that problem (2.1) is equivalent 
to the following one 



-r^u(x) - q(x)u(x) + Xu(x) - N(u(x)) =0, x G (0, l)\{a}. 
dx z 

u(0) = u(l) = 0; u'(0) = l, u'(a + 0) -u'(a-O) = /3u(a). 



(2.5) 



Atkinson (see [3], Ch.ll) proved that linear Sturm-Liouville problem with a potential 
as a function of bounded variation preserves the following properties of the Sturm- 
Liouville problem with a continuous potential: 

(a) all eigenvalues are real, simple, and form a monotone sequence increasing to 
infinity; 

(b) the sequence of the corresponding normalized eigenf unctions forms a complete 
orthogonal system in £2(0, 1). 

Also Atkinson have obtained asymptotic formulas for eigenvalues and corresponding 
eigenfunctions. 



In 32 ■ 37 the Sturm-Liouville problem with a potential as a function of bounded 



variation is studied using the generalized formulation of the problem. Particularly, 
in 32 35 36 it is shown that an eigenvalue A„ as a function of a potential q is 



analytical, monotone, and linear with respect to constants. 
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3 FD-method for solving Sturm-Liouville problem 
with potential including ^-function 

3.1 General description of the method 



Let us consider the following generalization of problem (2.5) 



£^u(x, t) - [ T q{x) + A(r)] u(x, t) - tN(u(x, t)) = 0, x £ (0, 1), x ^ a (3.1) 
u(0,t) = u(1,t) = 0; u'(0,t) = 1, u(a + 0, t) - u(a - 0, r) = 0, (3.2) 
<(a + 0,r) - <(a-0,r) = pu(a,r), Vr € [0, 1] . (3.3) 

The fact that for r = we can easily find an exact solution (A(0), Ui(x, 0), i = 1,2) 
to problem (3.1) - (3.3) and that for r = 1 the solution (A(l), Ui(x, 1), i — 1,2) to 



problem (3.1) - (3.3) coincides with the exact solution to problem ( |2.5| ), suggests an 
idea to write the solution to problem (3.1) - (3.3) in the form of series with respect 
to T 



u(x 1 t)=J2%V ,A(t)^It', Vz,tg[0,1]. (3.4) 
To apply the FD-mcthod's technique to the problem we also have to assume that 



(0 



d ^ d (%) i 

i=0 



dx 



—u(x,t) = }^— u(x)t 



dx 



i=0 



(3.5) 



for all t G [0, 1] and for almost all x £ (0, 1). 

Setting t = 1, we obtain the following representation 



V) 



A = A(l) = A , "0) = u(a;, 1) = ^ ( u(x), i = 1, 2 (3.6) 
j=o j=0 

provided that these series converge. Thus, we can represent the approximate solution 
to problem (2.5) as the pair of corresponding truncated series 



(J) 



i = l,2, 



(3.7) 



which is called the approximation of rank m. 

Supposing that the function u(x) £ C[0, 1] ( 3.6 1 satisfies boundary conditions ( 3.2 ) 

(i) 

we arrive at the conclusion that the functions u[x) are continuous on [0, 1] and satisfy 
the conditions 



u(0) =u(l) = 0, VieNU{0}, 



d ^u(x) 
dx 



d ii(x) 
dx 



x=0 



0, Vi G N. 



(3.8) 



a:=0 
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To meet requirement (3.3) we have to demand that 



d u{x) 



dx 



d u{x) 



dx 



x—a+0 



j3u{a), VieNU{0}. 



(3.9) 



x=a — 



Combining assumptions ( |3.4[ ), (3.5) together with equation (3.1) we make the 

conclusion that unknown pairs ( A (x)) , i € N U {0} can be found as the solutions 
to the following recurrence system of second-order differential equations 

d 2 (i) (0)(i) (i) 

u(x)+ A u(x) =F\x), x G (0, 1), x ± a, i e N (3.10) 

(*-p)(p) 



dx 2 

i-l 



F{x) = -2_^ A u(x) + q(x) u (x)+Ai_i IN; u(x), u(x ),..., u (x) J 



p=0 



d 2 (o), , (o)(o) 



(fx 2 



w(af)+ A u(a;) = 0, x £ (0, 1), x ^ a 



(3.H) 



supplemented with conditions (3.8), 
well-known Adomian polynomial (see 
N(u) 



3.9), where A k (iV; Uq, Wi, . . . , Vk) denotes the 



A k (N;v ,vi, ...,v k ) 



31 , [11]) of the order k for nonlinear function 
d k 



dt k 



,p=0 / 



t=0 



3.2 Basic problem 

(°) (o), , 

Let us consider the problem of finding A and u [x) in more detail. We call this 
problem the basic problem: 



d 2 (o), N (0)(0), 
dx 



2 u(x)+ A u(x) = 0, x € (0, 1), x ^ a, 



m(0)=u(1)=0, — 
dx 



x=0 



d u(x) 



dx 



d u[x) 



dx 



x—a+0 



/? u(a). 



(3.12) 
(3.13) 

(3.14) 



x=a—0 



(°), 



The unknown function u(x) satisfying problem (3.12), (3.13) can be represented 
in the following form 



(0) (0) (0) 

u(x) = u(x, A ) 



in (y\Wx\ x e [0,a] 

(VAM(l-i)) i6(a,l] 



(0) 

c sin 



(3.15) 
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where c e R. Using condition (3.14) and the fact that u(x) is continuous on [0, 1] we 

obtain the following system of transcendental equations for determination of A and 

(o) 
c : 



\ ,. /(0) (0) /. MU) \ 

sm[\l AaJ/V A = c sin(^V A (1 - a) J 



(°). 



(3.16) 



/(o) / A°) \ / A°) \ / A°) \ A°) 

V A cosfV A(l-a)J -cos(^V A aj = /3sin(^V A aJ/V A ■ (3.17) 



eg) 



(0) 



Eliminating unknown parameter c from system (3.16), (3.17) we arrive at the equa- 

(°) 

tion with respect to A 



/(0) 
A sin I 



t / Ao) \ / /(o) \ 

(V A J =-/8sinfV AaJsmfV A (1 - a) J (3.18) 



(0) 



and at the following expression for c 



(0) (0) (0) 
c = c ( A) = < 



(O) 

A a 



(0) 

A sin 



>s(/F) 



/ /(0) X 

(V A(l-a)J ^0, 
[V A(l-a)J =0. 



(3.19) 



Theorem 3.1. Suppose that a £ (0, 1) and k £ N, t/iera i/ie interval [it k ,tt (fc + 1) ) 
contains precisely one root of equation (3.18), provided that ft > 0. 

Proof. Suppose that /3 > and fix some positive integer k. 

First of all we consider the case of ft = 0. It is easy to check that in this case 
equation (3.18) possesses the countable set of solutions 



(0) „ „ 

A = 7r 2 fc 2 , VieN. 



(3.20) 



There are no other solutions of equation (3.18) except those presented in (3.20 We 
have that assertion of the theorem holds evidently. 

From now on we assume that ft > 0. For the given fixed k £ N there can be only 
two possibilities: 1) sin(7rfca) = 0; 2) sin(-7r£;a) ^ 0. 

(0) 

Let us consider case 1). We have that in this case the value A= 7r 2 fc 2 satisfies 
equation |3.18 1. Hence, to prove the theorem it remains only to show that the interval 
(ir 2 k 2 , 7r 2 (fc + l) 2 ) contains no other roots of equation (3.18). For this purpose we need 
to modify equation (3.18) like the following 

V 



f(y), y£ (7rfc,7r(fc + l)), 



(3.21) 



, (o) 
We do not take into account solutions with A < 0. 
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where 



f(y) 



sm(yaj sin^y(l — a) J 



sin 



w 



It is not hard to verify that 



f'(y) 



asin 2 (y(l — a)) + (1 — a) sin 2 (ya) 
sm 2 (y) 



(3.22) 



0. Equality sin(7rfca) = implies that sin(7r(fc + l)a) 7^ and we 
f(y) = +00. Taking into account the positiveness of f'(y) on 

1-0 

(nk, n(k + 1)) (see ( 3.22[ )) we arrive at the conclusion that (nk, n(k + 1)) — > (0, +00). 



This means that the graph of function z = — % can't intersect the graph of function 



z = f(y) on the interval (nk, n(k + 1)) (see Fig. [TJ and equation (3.181 has no roots 
on the interval (n 2 k 2 ,n 2 (k + l) 2 ). For case 1) the theorem is proved. 




a) 



Figure 1: The graphs of the functions 
z = f(x) (3.21 1 and z — —y/ ft, with a = 
1/3, f3 = 15. Equation (3.21[) possesses no 



roots on the interval (7rfc,7r(/c + 1)) with 
k = 3. 
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b) 



Figure 2: The graphs of the functions z = 



f(x) (3.21 ) and z = —y//3, with a = 1/5, 
/3 = 15. Here y* denotes the solution to 



equation (3.21 ) on the interval (irk, ir(k- 
1)) with k ■■ 



3. 



It is worth to emphasize that case 1) is possible if and only if a = t, where I € N, 
I < k. 

And now we switch our attention to case 2) sin(7rfca) ^ 0. It is easy to see that in 

(o) 



this case the value A = 7r 2 fc 2 do not satisfies equation (3.181. We intend to show that 
equation (3.18) possesses precisely one root on the interval (it 2 k 2 , 7r 2 (£; + l) 2 ). For this 



purpose we will use equation (3.211 again. As it was established earlier (see formula 
(3.22)) the function z = f(y) is strictly increasing on the interval (Trk,ir(k + 1)). 
Furthermore, it is easy to check that 



J 1 ™ f(v) 



-00, lim f(y) 

y-s-7r(fc+l)-0 



+00, if sin(7r(fc + l)a) ^ 0, 
0, if sin(7r(A: + l)a) = 



7 



and we have (7rfc,7r(fc + 1)) — > (— oo,0). The last fact provides that the graph of 
function z — f{y) necessarily intersects with the graph of function z = — -| on the 
interval (wk, n(k + 1)) in a single point y = y* , see Fig. [2] This point is a unique 

(0) 



solution to equation (3.21) on (Trk,ir(k + 1)) and A = (y*) satisfies equation (3.18) 



Hence, assertion of the theorem holds true, which was to be proved. 
The proof is completed. 



□ 



Proving Theorem |3.1| we have established an interesting property of the basic 
problem's solutions in regard to parameter a. We can formulate this property as the 
following remark. 



Remark 3.1. If a 

for all k e N. 



m (o) . . 

— , m, n £ N, m < n, then A = irkn satisfies equation (3.18) 
n ' ' 



Theorem 3.1 allows us to enumerate the roots of equation (3.18). We will denote 

(o) 



by X n the root of equation (3.18) lied in the interval \n 2 k 2 , 7r 2 (fc + l) 2 ). We also will 
use the notation (see (3.15), (3.19)) 



(0) , , (0). (0) (0) (0) (°) 
U n (X) =U{X,\ n ), C n =C(\ u ). 

Before passing to the next section we should make the following remark. 

(o) 

Remark 3.2. The value c n , n € N is always nonzero and its value squared can be 
represented in the following form 



/(o)y 



(o) 



( 



/ / (°) \ / IW) \ p f\l y "' \ 

[\l Kaj + cos[\j \ n aj H — sin^y X n aj 



1(0) 



(0) 



V 



(o) 



(3.23) 



4 Convergence of the FD-method 

In this section we will investigate the question of sufficient conditions providing the 



convergence of series (3.6). To obtain such conditions we will use the method of 



generating functions and the main part of this section is devoted to the derivation of 



the appropriate estimates for the terms of series (3.6) 



4.1 Equation for the generating function 



(o) 



Suppose that A n > is an arbitrary eigenvalue of the Sturm-Liouvillc problem (3.12 ) 



(3.14). The pairs A n , u n (x) for j = 1 . . .m can be found as the solutions to the second 

(4.1) 



order differential equations (see (3.10)) 
d 2 



dx 2 



U) (0) U) w 

u n (x)+ X n u n (x) =F n (x), x £ [0,1], i/a, 
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n [x) = - } X n u n {x) + q{x) u„ (x) + Aj-i[N;u n {x),u n (x), . . . , u n (x) 



p=0 



with boundary conditions ( see (3.8)) 



u„{0) =u n (l) = 0, — 

ax 



= 0, 



x=0 



and matching conditions ( see (3.9)) 



(4.2) 



d ul(x) 



dx 



d iil(x) 



dx 



x— a+0 



P u n (a). 



x=a—0 



(4.3) 



The general solution to the j-th equation of system (4.1) can be represented in 
the following form 



/ /(o) V 
sin I y\ n (x-€)} (j) 

-^-FniOdt, x€[0,a] 



U), , 
u n (x) 



(4.4) 



i / / (0) \ 
j- sin I Y X n (x - £)J 



Ci), s (i) . A/? 5 /! x\ 
u n (x) =c„ sin I y A„(l —x)\ 



U) 

Fn(Odt i6[o,l]. (4.5) 



For equation (4.1) to be solvable, the orthogonality condition 

l 

(i) (o) 

F n (0 «n (0« = 



(4.6) 



■ — . (i) 

have to be provided. Condition (4.6) allows us to find Ar 



3-1 



°r p) /" (°) ,^ (p) 



P =i 



X] A » / U n (£) U n 



(4.7) 



+ / °(0 ) (£) u n (£)df + / Aj-_i f JV; w°i (a;), u„ (a;), . . . , ( u„ ) (x)) u°„ , 



o 

where 



Mr, = / 

(o) y 

A„ o 



a . 2 

sin 



A/ (0) \ 
( y \ n x J 



Cnj / sin 2 ( y A„(l - x)J 



dx 
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/ (0) 

I ( fx sin(2yA„2; 
2 ~ 



(0) 



/ / (0) \ 

fio) 
4VA„ 



1(0) 

4\ A„ 



(o)\ 2 / a; 
' '"' 2 



sin 



1 / s,:l 

2 A " 2 V A 



in(^2YA„aJ^ (c„J ^ ^ sin(2y A„(l - a)J 



/(o) 
2\/ A r , 



It is easy to verify that system (3.16), (3.17) implies the equality 



( / (0) \ ( / (0) \ i( / (0) \ 

(0) ^ 2 sin( 2y A„(l - a)J sin(2yA„(a)J /3sinfyA„a) 



(c°„) 



1(0) 

2 V X n 



(o) /(o) 



(0) N 2 
A n 



using which we can obtain the following representation for M n : 



Mr, 



(0) 

An / , 0) , /(o)\ 2 



\ - „ /, / (°) \ /W\ A \ 

a I A„ + ( c„ (l-a)+/3sin^ VA„a / A„ . (4 



On the other hand, formula (3.23) can be simplified like the following 



/(0)\2 _ 

V n ) ~ (o) 

An 



H i= cos^y A„a j sm^y A„a J + — sin A„a J 



'(0) 



(0) 



P 2 . ,, 

An 



1(0) 



(4.9) 



Combining formulas (4.8) and (4.9) we arrive at the following expression for M n 

f(o) 



M„ 



1 



1 



(o) N (0) 



/3 sin (2 V A„aJ / V A„ + (3 2 sin 2 V X n aj / A 

9 A A ^ \ , (o)N 

+/3sin 2 (^\/A„aJ/ A,, 



(l-a)+ (4.10) 



or 



For the function u n (x) (|4.4|, (4.5) to be continuous in the interval [0,1] and 



CO . 



satisfy condition (4.3), the parameter c„ in formula (4.5) should satisfy the following 
equations 

! f Jm \ 

(j) . fJ^K, ^^ f sm [ \ ^n{a - 0) U) 



u > ■ A/wi \"\ 
c„ sml y A„(l -a) 



(4.11) 
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, <v<> \ r cosl V A„(a-£) o) /• sinl V A„(a - 0) o) 
£ cos(V A„(l-a)J = - / ^ Frffidt-p I V ^ ^ F^)de- 



'(0) 



'(o) 



(0) 
An 



(4.12) 

Formulas (|4. 11 [) and (|4.12[) will be used in the software algorithm of the FD- 



method, as it is described in section [5j However, to obtain the convenient estimate 



for |c„| we have to square both sides of equalities (4.11) and (4.121 and then sum 
them up: 



\ c n\ ^ I A, 



(i) 



1(0), 



FrlO di+{P/\j\ n ) J FrlO 




U) 



1/2 



< 



(0) 



1 + (1 + /3/V A ri ) 2 } 



< 



FrlO 



' /(0) 

l + (l + /3/VA„) 5 



(o) 

Frlx) . (4.13) 

0,1 



Using estimate ( |4. 13 1 and formulas (4.4), (4.5) we can easily estimate the value of 



\un\\ = max |i<ia;)| like the following 

1 llo ° a;G[0,l]' 1 



(j) (0) 



(4.14) 



rJ- 1 



< a 



2J A ™ llr™lloo + IMIo,i II u n l| 00 + ^-i( iV ;|r™lloo'---'ll u « I 

p=0 



where ||q||o,i = / \q{x)\dx and 
o 



a = a(n) — 



' /(0) \ /(0) 

l + (l + /?/VA„) 2 + l /VA n , 



(4.15) 



_ oo (0) 

and iV(w) = ^2 \ a p\ uP ■ Similarly to that we can estimate the value of |A n | (see 



formulas (4.7), (3.15 1) 



Aj < 6 



rJ- 1 



(j r P) Ml(P)|l 0'-l)ll „ (Zr l|(°)|l I|CJ-1)| 



Y^l \ III (P) II _L_ 4 (AT 

ZJ A ™ IIKIloo + Nlo.l II U ™ lloo+^J- 1 ^ 
P=l 



(4.16) 



where (see (3.23)) 



f ■' 1 " ' , / — 1 1 " ' i 1 1 1 1 • 
b{n) = maxjl, \j \ n /M n , \j \ n \/c n /M n j >A„ \\uXx)\\jM n . (4.17) 



(0) 
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Using notation 



b I, W 

tin 



/iiU „ l!x and N = \X n \/a j -\ j G N U {0} (4.18) 
we can rewrite inequalities (4.14), (4.16) for j = 2, 3, . . . in the following form: 

3-1 



where 



Vj < ^2v<3-pv p + ||g||o,iUj-i +A 7 --i(iVi;0, Uj_i), 

3-1 

Mi < y^Vj-p v p + lkllo,i«i-i + A_,-_i(JVi;0, ui,... .w^-i), 
P =i 

iv 1 («) = iv(||^|| oo + «). 



(4.19) 



To obtain similar estimates for V\ and \X\ let us consider inequalities (4.14), (4.16) for 
j = 1 in more details. Thus, 



|(i). 



< a 







(|A n | - 


HMIo.l) 



lkl|o,l||w„|L+^(h nL) 



(4.20) 
(4.21) 



Using the fact that N(0) — and inequality (4.20) we obtain the following estimate 
of ui : 



(4.22) 



«i < (jii + ||?|lo,i) v o + bN(\\ul\\ 

= (Mi + lkll 0l i)^ + M^(l|SilL)-^(°)) 

= (/ i l + lklloa) Wo + 6 ^'( |l"«lloo)ll'"™lloo - 

< (mi + IMIo.i) v o + N'(\\un\\ 00 ^v a =(jn + ||g|| 0)1 ) Vo + N[(0)vq. 

Similarly to (4.22) we can get the following inequality for fi^: 

Mi < lkllo,i + iV((0)^ O - 

Let us consider two sequences of positive real numbers {%}°1 and {^1°% de- 
fined through the following recurrence formulas 

(4.24) 



(4.23) 



Mi = IMIo.i v o + N[(0)v Q =vi- favo, 
vi = (mi + IMIoi)«o + ^(0)«o, 



(4.25) 
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j-1 



= 5Z' i i-p w P + Wl\\o,i v j-i + Aj-i(Ni,Q,vi, ...,Vj-{) = Vj - /j.jVo, (4.26) 



3-1 



p=0 



(4.27) 



where (see (3.23), (4.17) and (4.18 1) 



l(0)l 



< max-^ 



= = v (n) = bin^lunW^ < 
/iui . /(()) , , f ,/(0) 



r / w / w i r / w < / ^ i 

jl, Y X„/M n , Y A n V c n/Mn|max|l/y A„, y/c n /y A„j 



(4.28) 



= maxjl, \/c^| maxjl, y^, M„/y A„|/M n . 
Taking into account that the inequality < 1 implies the estimate M„ < 1, we get 



from (4.28) the following estimate for vq : 



v < maxjl, d, vc„M„/Y A„ j/M n . 



(4.29) 



It is easy to see that Vj < Uj an d /i 7 - < ^ Vj € N. 

Eliminating from system (4.24) - (4.27), we obtain recurrence formulas 



«i = (l + wo) (lkllo,i«o + ^i(0)«o) , 



(4.30) 



4-1 



(l + »o) (||g||o,i^-i+^j-i(JVi.O,« 1 ,...,l;j-i)) , i e N\{1}. 

(4.31) 



Denoting by /(z) the series (generating function) 



/(*) = £' 



(4.32) 



j'=i 



and using equalities (4.30), (4.31) together with iVi(/ (z)) = ^ z J Aj ( JVi ; 0, iJi , . . . 

3=0 

we arrive at the equation with respect to f(z) : 

f(z) = [/(*)] 2 + z(l + v ) [||g||o,i(/(«) + v ) + Ni{f{z)) + N[(0)v ~ Ni(0)] ■ (4.33) 
4.2 Convergence result for the linear case (N(u) = 0) 



Let us first consider the case when N(u) = 0. Equation for generating function (4.33) 
can be rewritten in the form 



f(z) = (f(z)) 2 + z(l + v )\\q\\oAf(z)+v ) 
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or, which is more convenience, 

(/(z)) 2 - [1 - (l + v )\\q\\o A z\f(z)+v (l + v )\\q\\ OA z = 0. 
We have the quadratic equation with respect to f(z) with the roots 

(l-(l + 1;o)||g|| 01 z)±V^ 

Adz) = ^—^ > 

where 

D = (w 1 - w 2 z)(l/w 1 - w 2 z), 

w 1 = 1 + 2v + 2y%o (1 + v ), u; 2 = (l + wo) Ikllo.i- 
The solution which represents the generating function is 

, i-(i + «o)||g|lo,i»--/D 

/ ( Z ) = n ' ■ 



(4.34) 



(4.35) 



It is obvious that the right-hand side of equality (4.35) can be expanded as a power 
series in z, Vz 6 [0,R] : 



v 1 (l+^o)IMI ,i 1 
/(2)= 2 2 - Z ~2 



w 2 ^ (2p-3)H i-p p ' 
2^ ^ (2p)H 1 2 



x (4.36) 



1 w 2 y/m ^(2p-3)!! -i+p ' 

7. z — S 77. \ . i 1Ui M),Z f 

2 ^ ( 2p )!! i 2 



i (i + «o)||g|| ,i 



l l 

-2 + 2 



(2j - 3)!! 
(2i) 



j'-i 



irH+v)-£ 



(2p-3)!! (2j-2p-3)!! 
^ (2p)!! (2j-2p)!! ' 



-i+2p 



w 2 ((tui 



i=2 
where 



(2j - 3)!! 
(2j)» 



(l+< 2j ) 



-2+ 



(2j)!! U(2p-3)!!(2j-2p-3)! 



HI E 



,- 2 P 



(2j-3)!!^ (2p)!! (2j-2p)M 1 



i? = i?(n) 



WiW 2 



l + 2v Q -2y / v (l + v Q ) 

(i+«o) Ikllo.i 



t; = «o(n), 



(4.37) 



dof 



(2fe)!! = 2-4-...-2fc, (2fc + !)!! = !• 3-... • 2fc + 1, (-!)!! = ! 



Taking into account that series representation (4.36) is valid for z = R, we arrive 

ients 

(K 



at the following inequality for coefficients of generating function (4.32) 

-l 



< Rvi 



R<1 
~ 4 
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[)<K ^- 2 (2i)!! 



1 + 1% 



2j 



(2j)!! ^(2p-3)!!(2i-2p-3)H 



(2j-3)!!§ (2p)!l (2j-2p)H ^ 



-2p 



< 



(2j-3)!! . 

Using the Stirling's formula, we can estimate «j like the following 
(2j-l)!! _ (2j)! _ (2j)! 



(4.38) 



3 2(2j-1)(2j)!! 2(2j - 1) ((2j)!!) 2 2 2 J+ 1 (2j - 1) (f-f 
2y/¥] i^f 3 eW 1 



< 



(4.39) 



2^(2i-l)(^(i) J j 2(2, - 1)^7 " (2, - 1)^7 
Returning to notation ( 4.18| , we arrive at the estimates (see (4.17) and (4.29)) 



I I 
\u n \ 



a J Vj a 3 _ 

< — Vj < mini 

b b 



{l, M„/V A n , M„/(y S vC) Ja.-Ti (4.40) 



and 

U) 

| A„ | = a? fj,j < a 3 < 



VjO? 1 



< 



Q!jM„ 



( M„ + maxjl, M n /]/Zy\R 



where 



i+yi+(i +( 8/yA n ) 2 



(4.41) 



(4.42) 



Inequalities (4.40), (4.41) imply the error estimates 



^ „())„ f /(°) / /w , — \i 

^ H P"L <Cn,mmin[l/M n ,l/YA n ,l/(V^n^)} 



tin— ti 



m+1 



V" l W l 

|A n — A ri | < 7 ^ | A n | < 

j=m-\-l 



( M n + maxjl, 

av\A«}) 



1? 



provided that 
where 



r n < 1, 



1 - 



(4.43) 

(4.44) 
(4.45) 



i+V i+(i+/3/y A„)2 

(0) 

A„-R 
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Theorem 4.1. The FD-method described by formulas (3.7 1, (3.15), (3.18), (3.191, 



(4.4), (4.5), (4.7), (4.11), (4.12) converges superexponentially to the exact solution of 



problem (2.1 1 with N(u) = provided that condition (4.44) holds. The error estimates 



for the convergent FD-method are described by formulas (4.28), (4.42), (4.43), (4.45) 



with R determined by formula (4.37) and aj determined by formula (4.38) 



3.1 



(0) 



it follows that A n > ir 2 n 2 . Using this fact, it is not 
(4.10) tends to 1/2 as n — > +oo andc n (3.23) tends to 1 as 



Remark 4.1. From Theorem 
hard to verify that M n 

Thus, we can conclude that for fixed a e (0,1) and j3 > there exists a 



-oo. 



positive integer no such that for all n > uq condition (14441 holds true 



4.3 Convergence result for the nonlinear case 

In this subsection we derive the conditions providing the convergence of the FD- 
method for the general case when N(u) ^ 0. 

Equality (4.33) holds for any z in the convergence interval of the power series 



(4.32) with the coefficients Vj determined via recursive formulas (4.28), (4.30), (4.31 ) 



By R we denote the radius of convergence of this series. We intend to prove that R 
is nonzero (positive). For this purpose let us consider the inverse mapping z = f . 
From (4.33) we obtain 



*(/) 



f-f 



(1 + u ) [Nlo,i(/ + «o) + Nx(f) + N[ (0) v Q - m (0)] 



(4.46) 



Since z(f) is holomorphic in some interval containing the point fa = and z (0) = 0, 
we can directly calculate the derivative z' (0) 



Z '(0) = 



lim 

/->o 



*(/)-*(0) 



f-o 



(4.47) 



lim — 

/->o (1 



1-/ 



«o)[lkllo,i(/ 
1 



v (l + v )[\\q\\ 0!l +Ni(Q)] 



v Q )+N 1 (f) + N[(0)v 
> 0. 



Ni (0)] 



Inequality (4.47) implies that there exists an inverse function / {z) which is holo- 
morphic in some interval (—R,R) [9) p. 87]). Now let us prove that series (4.32) also 
converges at the endpoints of the interval. For this purpose it is enough to consider 



only right endpoint z = R. Conversely, suppose that series (4.32) diverges at the point 
z = R, that is, 

lim f(z) = +oo. 



However, taking into account that equality (4.33) holds for all z in (—R,R), we get 
the following contradiction 



1= lim /(*)+ (4.48) 

z~ >_R— \ 



(1 + i7o)z[||g|lo,i(/(z) + go) + mifjz)) + N[ (0) v - Nx (0)] 



+oo. (4.49) 
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This contradiction implies the inequality f(R) < +00. Thus, for some positive con- 
stants c and s the following inequality holds 



R>V4 < — 

i 1 ' 



where the constant R = R(n) is uniquely determined by ||g||o,i) vo = v o( n ) an d Ni(u) . 



Returning to notation (4.18) we arrive at inequalities (4.401, (4.41) and error 



estimates (4.43), as it was in the linear case. However, in the nonlinear case, parameter 



R denotes the convergenc e rad ius to the power series (4.32) satisfying equality (4.33) 
and ctj = c/j 1+£ . Remark 4.1 is consistent with the nonlinear case too. Taking it into 



account we arrive at the following theorem. 

Theorem 4.2. For any a G (0, 1), @ > 0, q(x) € Li(0, 1) and N(u) £ C°°(R) there 



exists a positive integer no, such that Vn > rio condition (4.44) holds and the FD 
method described by formulas dO), (ET51), p38|, pT9|, (44L (gjL dO), OT) 



(4.12) converges superexponentially to the exact solution of problem (2.1). The er- 



ror estimates for the convergent FD-method are described by formulas (4.43), (4.42) 
where R denotes 



(4.45), (4.281 



the 



convergence radius to the power series (4.32) 
satisfying equality (4.33) and otj — c/j 1+£ for some positive constants c and e. 



5 Algorithm of the FD-method: software imple- 
mentation 

In this section we discuss the question of possible software implementation for the 
proposed FD-method. 



Generally speaking, to apply FD-method to problem (2.1 ) we should perform two 



main steps: the first one is to solve the basic problem (3.12) - (3.14) and the second 



0') 



(j) 

one is to calculate a certain number of corrections A„, u n (x), sufficient to achieve 
the required accuracy. The second step includes the repeated application of formulas 



(4.4), (4.5), (4.7), (4.11 ), (4. 12 ), which contain the operator of integration. Hence, the 



second step can be executed analytically (using analytical integration) only for the 
rare special cases, and in practice we almost always have to use numerical integration 
instead. Using numerical integration, however, we approximate cigenfunction on the 
quadrature mesh only, whereas analytical integration yields us the uniform approxi- 
mation on the interval (0, 1). In the algorithm described below we approximate the 



integral operators using Sine quadratures and Stenger's formula (see 29 



K 



(6-0) 



f(x)dxKh s Sk-vf ( ~r — r~ . 

n ' p £^ K PJ \ l + ePhs J ( e - phs /2 + ePhs /2y 



(5.1) 



where 



a + be h ° k 1 
z k= l + e h sk > k = -K...,K, S k = - 



sin (irt) 



nt 



dt, k = -2K ...2K, 



3 The function f(x) is needs to be sufficiently smooth on (a, b), see [29]. 
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and h s = for some < d < ir, [i > o|*| 

Also, we make the assumption that function q(x) is analytical on the intervals 
(0, a), (a, 1) and might has the integrable singularities at the points 0, a, 1 only. 
However, the proposed algorithm can be easily generalized on the case when q(x) has 
any finite number of the integrable singularities on (0, 1). In the reminder part of the 
section we will use the notation 



1 + e h *j 



Z-2. 



a + e 



h.j 



1 



oh,j 



(5.2) 



Mi,. 



a 



(e-jh e /2 _|_ e jh e /2y 



M2, : 



( e -jh./2 + e i h -/ 2 Y 



(«£>(*)) ' 



j = -K,...,K. 



Evidently, the highest FD-method's accuracy which can be achieved using Algorithm 
[l] presented below is restricted by two factors: the accuracy of quadrature formulas 
and the accuracy with which the basic problem can be solvecj^] To emphasize this 



fact, we will denote by A^ p ' the value of A n perturbed with this two factors. For the 
same reasons we will use the notation U™'^ instead of u n {zi t j), notation instead 



(p) 



4 For discussion on the optimal choice of value for h s , and parameters d and see [29] , [30| . 

5 Here we intentionally do not consider yet another factor: the machine precision, we assume that 
it is high enough (arbitrary). The arbitrary precision in the computational software can be provided 
via the libraries for multiple-precision floating-point computations, like MPFR |§J. 
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of A„ and notation U i j instead of u n (Zi t j) 



Algorithm 1: FD-method's algorithm for solving problem (2.1 1 



input : The real numbers a g (0, 1), (3 > 0, e > 0, functions q(x) <G £i(0, 1), 
N(u) E C*°°(R), N(0) =0; 

the number n > of eigensolution to be approximated; 
the rank r of the FD-mcthod; 



output: The value such that |A„ — A^|< e and the array of values U™'^ , 
such that \u n (ztj) — f/"j r | < e, i = 1, 2, j = —X. . . . , K; 

begin 



SolveBasicProblem(rt); // Solve the basic problem 

FindParametersForQuadratureFormulas; 

for k := 1 to r do 

FindNextCorrectionForEigenvalue(fc);// Find An 
FindParameterC(fc);// Find 



FindNextCorrectionForEigenfunction(fc);// Find U i 



,(fe) 



end 



K ■■= E Al p) ; 

for j := — K to if do 

,(p). 



p=0 p=0 



end 
end 



Let us consider Algorithm [T] in more detail. 

Procedure SolveBasicProblem(n) calculates the approximation for solution to 



(0) 



<?>_ MO) _ 



the basic problem ( |3.12[ ) - ( |3.14| . If na <E N then A„= A^ ' = ir 2 n 2 and c„= C, 

(o) 



cos(irn)/\J = (— l) n+1 / \/ A^; if not, then the approximation A„ ^ for A„ can 



(0) 



be found from equation (|3.18|) via Newton's method and approximation Cn " 1 for the 

r n,(0) 



parameter c„ can be found using formula (3.19). After that, the approximation U, 



(0), 



for u n (zij) can be calculated using formulas (3.15), (3.19) as it is described in the 
procedure. 

Procedure FindParametersForQuadratureFormulas calculates the values of 



parameters h s and K, such that when quadrature formula (5.1 ) is applied to integrals 



(4.4), (4.5), (4.7), (4.11 ), (4.12 ) it provides the required accuracy e. Actually, the pa- 



rameters h s and K can be found through the a posteriori error analysis of quadrature 
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formula (5.1| applied to the integral J (q(x)+ u^x) + N($l(x)^dx (see H)- 



Procedure SolveBasicProblem(n) 



input : The number n > of eigensolution to be approximated; 

(0) l (0) (0)1 

output: The value A„ ; such that A„ — A„ ; < e and the array of values 
U?f\ such that \u ] n (z h3 ) - Ulf ] \<e,i = l,2,j = -K...,K; 



begin 



if na € N then 



Af := nn; CT := (-1)«+7A 



-,(()) 



n+l/\(°). 



else 



A^ 0) := nn; A : = 0; 
while |Al 0) - A | > s do 



// Using the Newton's method 

s do 

A := Al 0) ; Ai 0) := A - // Here /(A) = . n ™™ „ + f 



end 

r-(o) ._ 



i(AWq) 



A™ sin(Al '(l-a))' 



end 



(0) 



// Compute approximations for u n using formulas (3.151 

for j := —K to K do 

I U?f := sin^z^/A^; U?f := sin(A^(l - z 2j )); 



end 

A< 0) := S<Z r(A< 0) ); 



end 



Procedure FindNextCorrectionForEigenvalue(fc) calculates the approxima- 



tion An for A„ using formula (|4.7|). The integral in the right side of this formula is 



(k) 



(k) 



approximated using tanh rule, see 22 . Also, this procedure computes the approxi- 



mation Fij for F n (zij) (see formula (4.1)) 



Procedure ComputeParameterC(fc) computes the approximation C„ for pa- 
rameter c„ using formula (4.11) if an € N or formula (4.12) if an f N. In the case 



when an € N formula (4.12) was used, though, in a slightly simplified form. Taking 



(0) 

into account that in this case A n = in, this formula was rewritten as following 



c„= (-1) 



cos(7r<) (j) rsin(7r<j (j) 



Kit 



II 







/--R- (5.3) 

(nn) z 



Eventually, procedure FindNextCorrectionForEigenfunction(fc) calculates the 



approximation £/™j^ for it^(zjj) using formulas (4.4) and (4.5 ) together with Stenger's 



(fc) 
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formula (5.1 1. 



Procedure FindNcxtCorrectionForEigenvalue(fc) 



input : An ^ for p = 0, 1, . . . , k — 1; 

C/"j (p) for i = 1, 2, j = -K , . . K, p = 0, . . . , k - 1; 



output: A! F itj mF^Zij), for i = 1,2, j = -K ...K; 
begin 

for j := —K to K do 

fc-i 

*U := - 

P =i 
k-i 



(fc) 



If' 



2J 



end 

for j 



E A^Ulf+ q (z 2J )U^ k - 1) +A k ^(N-, Ul 
p=i 



(fc) 



(0) 



>(fc-l)<i 



2J 



A„ :=0;// Compute A„ (see (4.7)) via the tanh rule, see [22] 
-if to if do 



I a(*0 \( k ) , rr™.(0) p , rrn>M r ,, 

| A n .— l\ n + U 1 j rij/Jij + U 2 j ^2,j^2,ji 

end 

A (fe) ._ u a(*)/ . 

for j := —if to if do 

I p p AWr^C 1 ). p p a (fc)rr n '(°). 

end 



end 



Procedure ComputeParameterC(/c) 



input : A^ 0) ; F hJ for i = 1,2, j = -K ... AT; 



(A;) ■ 1 i 

output: Cn «c n , see formulas (4.111 and (4.12) 
begin 



A 



AW; Ck k) :=0; 



if na € N then 



// ...then A = 7r?i and we use formula (5.3) 
for j := —K to K do 

r (k) ._ 
<-^n ■ — 

d fe) - (cos(Azij) - ^sin(Azi J )/A) Fijfiij - cos(Az 2 , J )F 2 j^ 2 j; 
end 

Ci fe) := (-l) n h s C ( n k) /A; 



else 



// ...else we use formula (4.11) 
for j := —K to K do 

| C ( n k) := C$7 + sin(A(a - zij^Fijfj.ij + sin(A(a - z 2 j))F 2yj ^ 2j ] 
end 

Ci fe) :=/M?i fc V(Asm(A(l-a))); 



end 



end 
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Procedure FindNextCorrectionForEigenfunction(fc) 



input : Al 0) ; ci k} ; F itj for i = 1, 2, j = -K . . . K; 



n,(k) ( k ) 



output: U" k 
begin 

A:= v^; 



[Zi,j) for i = 1,2, j = -K. . . K, see jOJ and ((45 1 



9 
10 



12 



for p := —K to X do 

Zi := 0; Z 2 := 0; 
for j := —If to K do 

z i,p (fc) 

Tx :=li + cos(Azi t j)Fijfj,ij5p-j; // Ii« J cos(A£) F„ 

o 



Z 2 := Z 2 



o 

zi.p (fc) 

o 



end 

TT n,(k) 
u l,p 

end 

for p := —K to if do 
Zj := 0; 2 2 := 0; 
for j := —K to K do 



ft. s ^sin(Azi. p )Ii — cos(Azi !P )Z2^ /A;// See formula (4.4) 



l (fc) 

X 1 :=l 1 +cos[Az2,j)F 2i jH2,jSj- p ; II I 1 « J cos(A£) F„ 



2 2,p 
1 



1 (fc) 

Z 2 :=Z 2 +sin(Az 2j )F 2j p 2j <5 J _ p ; // Z 2 « J sin(A£) F„ (£R 

-22, p 



end 



17. 



-n,(fe) 



2,P 

c 



end 
end 



'n ] sin(A(l - z 2 . p )) - h s (sin(Az 2:P )Zi - cos(Az 2 , p )Z 2 ) /A 
formula p~5| 



;// See 



Using the proposed algorithm and the MPFR library |6j, a C++ application was 
developed for solving problems of type ( 2.1 ). The numerical example presented in the 



next section was prepared via this application. 



6 Numerical example 



Let us consider the Sturm-Liouville problem (2.1) with 

a = \, 13 = 2, q(x) = H — . 1 H — . 1 H — , 1 , (6.1) 
2 y/\Q.7-x\ v/|0.1-x| v/|0.3 -a;| A /|0.4-2a;| 

AT(u) = u 9 . 

For the error control we will use the following functionals 
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d u n (£) 



d£+ (6.2) 



1 -77 h P u„(a) + 



7T7 771 771 

A„(a) = 3^ u„(q + 0) - — u n [a - 0) — P u n (a). (6.3) 

To tackle problem (2.1) with (6.1) the algorithm described in section [5] has been 
slightly modified. The intervals (0, a) and (a, 1) have been split into subintervals 
(0, 0.1), (0.1, 0.2), (0.2, 0.3), (0.3, a) (a, 0.7), (0.7, 1) with respect to the singular points 
of function q(x). After that all integrals over the interval (0,1) were substituted by 
the sum of integrals over each subinterval. The approximations for the ten eigenvalues 
and eigenfunctions were successfully found and presented in Tab. [I] and Fig. [3]-[5j 
The graphs presented on Fig. |6j[7]confirm the exponential nature of the FD-method's 
convergence rate. 



Table 1: Example 1. The results of computations. 



n 


m 


Ttt 

An 


|| WnOzOlloo 


(m) 

1 A» | 


771 

r n 


A„(a) 


1 


10 


23.437363200234028176652 


1.5e-ll 


7.6e-10 


2.5e-ll 


2.4e-26 


2 


10 


50.879953432153777724296 


7.5e-12 


2.4e-10 


2.0c-ll 


5.6-27 


3 


10 


102.294039773949565868154 


3.0e-13 


5.7e-ll 


1.9e-13 


9.7e-27 


4 


10 


167.932111361326104363494 


2.2e-16 


1.2e-13 


l.le-16 


1.9e-27 


5 


10 


261.703789042290324125067 


2.5e-17 


6.8e-16 


3.3e-17 


8.3e-27 


6 


10 


365.290665054662412777331 


1.5e-17 


9.5e-16 


3.7e-18 


5.3e-27 


7 


10 


497.311217072847814939907 


6.3e-19 


1.2e-16 


l.le-19 


7.7e-27 


8 


10 


642.305601325675356973240 


6.7e-19 


4.6e-17 


1.3e-19 


l.le-26 


9 


10 


813.233561353244869046018 


7.3e-21 


1.7e-17 


9.1e-22 


8.5e-27 


10 


10 


995.761252385458344653891 


3.7e-21 


2.9e-18 


4.4e-22 


1.8e-28 



Also, to illustrate the fact that the method's convergence rate increases along with 
the number of approximating eigenpair we have computed linear approximations to 

/ ( m ) A fS m \\ 

the functions y u . n = y u , n (m) = ln(J| u n (x)\\ j , y x , n = yx, n {m) = ln( |A„|J, y r<n = 
Vr. n {m) = ln(r") : 

V\,n = a x . n m + b x , n ~ y\,n(™>)> 

These approximations have been obtained via the method of least squares using the 
values of corresponding functions at the points m — 0, 1, ... , 10. The parameters of 
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linear approximations together with approximation errors e. u ,ni e \,m e r,n (see (6.4)) 
are presented in Tab. [2] 

e«,„= max \a un m + b un -y un (m)\, 

rn— 0,1,. ..,10 

e A ,„ = max lawm + bx, n - yx.n(m)\, (6.4) 

m— 0,1, ...,10 

e r „ = max |o r , n m + b r n - t/ r n (m)|. 

m=0,l,...,10' 

From table [2] it follows that the coefficients a u .„, ax, n > a r,n decrease almost monotoni- 
cally as parameter n increases. Taking into account the comparatively small amounts 
of approximation errors, we can conclude that the convergence rate of the method 
has an exponential nature and it does increase as n tends to infinity. 



Tabic 2: Example 1. The parameters of linear approximations. 



n 








d\,n 


bx,n 


e\,n 




b r n 




1 


-2.3 


-1.7 


2.0 


-2.4 


2.8 


2.2 


-2.4 


-2.0 


2.1 


2 


-2.3 


-2.0 


0.6 


-2.3 


1.7 


5.7 


-2.3 


-2.8 


1.9 


3 


-2.6 


-2.4 


0.8 


-2.8 


3.5 


1.8 


-2.7 


-3.0 


1.3 


4 


-3.3 


-2.5 


0.9 


-3.5 


4.1 


1.8 


-3.4 


-3.1 


0.9 


5 


-3.5 


-2.4 


0.3 


-3.8 


4.9 


2.1 


-3.4 


-4.0 


0.8 


6 


-3.6 


-2.9 


0.6 


-4.0 


4.7 


4.8 


-3.6 


-4.1 


0.8 


7 


-3.9 


-2.7 


0.5 


-4.1 


5.0 


2.4 


-4.0 


-4.0 


0.4 


8 


-3.9 


-3.1 


0.6 


-4.2 


5.5 


2.2 


-3.9 


-5.0 


1.1 


9 


-4.3 


-3.3 


0.6 


-4.5 


5.2 


1.6 


-4.4 


-4.7 


0.9 


10 


-4.4 


-3.0 


0.7 


-4.7 


5.2 


2.5 


-4.4 


-4.7 


1.0 



On the other hand, the values of parameters Vo(n), (see (4.28)), R(n) (the con- 
vergence radius for generating function (4.32), ( 4.33[ )) and r n (see ( |4.44[ )) calculated 
directly for the case of problem (2.1 ), (6.1 ) indicate that the FD-method for eigenpairs 
with numbers n = 1, 2, . . . , 10 have to be divergent (r„ > 1), see Tab. [3], whereas, in 
fact, it is convergent. This means that the convergence conditions stated in Theorems 
|4.1[ |4.2| are essentially overestimated. 



Table 3: Example 1. The values of parameters vo(n), R(n) and r. 



n 


v (n) 


R(n) 


r n 


n 


Vo(n) 


R{n) 


r n 


1 


1.8 


0.41e-2 


189.9 


6 


2.0 


0.34e-2 


39.1 


2 


2.0 


0.34e-2 


125.1 


7 


2.0 


0.34e-2 


33.0 


3 


2.0 


0.34e-2 


76.5 


8 


2.0 


0.34e-2 


29.1 


4 


2.0 


0.34e-2 


59.6 


9 


2.0 


0.34e-2 


25.6 


5 


2.0 


0.34e-2 


46.3 


10 


2.0 


0.34e-2 


23.2 
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Figure 5: Example 1. 
(figure b)) 



FD-method. The graphs of 113 (x) (figure a)) and 
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Figure 7: Example 1. The graphs of the 

/ (™) A 

function y\. n (m) = In I A„ J for n = 
0,1,- --,4. 

7 Conclusions 

Summarizing the theoretical and practical results presented in the paper we can con- 
clude that the Sturm-Liouville problems with potential including <5-function can be 
successfully treated with the FD-approach. The authors of the paper are unaware 
of the software packages which can solve the problems of such type. However, the 
presence of (^-function in the potential do not introduce significant changes to the FD- 
method's algorithm in comparison with that for the classic Stourm-Liouville problems. 
The algorithm of the FD-method can be easily modified for the case of potentials with 
a finite number of integrable singularities, as it was described and discussed in sections 
[5] and [6] The numerical example presented in section [6] confirms the predictions of 
Theorems |4.1| and |4.2| about the exponential nature of the FD-method's convergence. 

However, as it was mentioned in section [6j the convergence conditions stated in 
theorems |4.1[ |4.2| are essentially overestimated and the search for the more subtle 
conditions is still a pressing issue. 

References 

[1] S. Albeverio, F. Gesztesy, R. H0egh-Krohn, and H. Holden. Solvable models in quantum 
mechanics. AMS Chelsea Publishing, Providence, RI, second edition, 2005. With an 
appendix by Pavel Exner. 

[2] S. Albeverio and P. Kurasov. Singular perturbations of differential operators, volume 
271 of London Mathematical Society Lecture Note Series. Cambridge University Press, 
Cambridge, 2000. Solvable Schrodinger type operators. 

[3] F. Atkinson. Discrete and continuous boundary problems. New York, 1964. 



27 



B.I. Bandyrskii, LP. Gavrilyuk, I.I. Lazurchak, and V.L. Makarov. Functional-discrete 
method (FD-method) for matrix Sturm-Liouville problems. Computational Methods in 
Applied Mathematics, 5(4):l-25, 2005. 

Ronald N. Bracewell. The Fourier transform and its applications. McGraw-Hill Series in 
Electrical Engineering. Circuits and Systems. McGraw-Hill Book Co., New York, third 
edition, 1986. 

Laurent Fousse, Guillaume Hanrot, Vincent Lefevre, Patrick Pelissier, and Paul Zim- 
mermann. MPFR: a multiple-precision binary floating-point library with correct round- 
ing. ACM Trans. Math. Software, 33(2):Art. 13, 15, 2007. 

LP. Gavrilyuk, A.V. Klimenko, V.L. Makarov, and N.O. Rossokhata. Exponentially 
convergent algorithm for nonlinear eigenvalue problems. IMA Journal of Numerical 
Analysts, 27:818-838, 2007. 

R.G. Gordon. New method for constructing wavefunctions for bound states and scat- 
tering. Journal of Chemical Physics, 51(l):14-25, 1969. 

Einar Hille. Analytic function theory. Ginn and Co., Boston, 1959. 308pp. 

Liviu Gr. Ixaru. Metode numerice pentru ecuatii diferne^iale cu aplicafa. Editura 
Academiei Republicii Socialiste Romania, Bucharest, 1979. With an English summary. 

Abbaoui K., Cherruault Y., and Seng V. Practical formulae for the calculus of multi- 
variable adomian polynomials. Math. Comput. Modelling, 22(l):89-93, 1995. 

N. Kryloff and N. Bogolioubov. Sopra il metodo del coefficient constati (metodo del 
tronconi) per l'integrazione approssimate delle equazioni differenziali delle fisica math- 
ematica. Bolletino della Unione Mathematica Italiana, 7(2):72-76, 1926. 

V.L. Makarov. About functional-discrete method of arbitrary accuracy order for solving 
sturm-liouville problem with piecewise smooth coefficients. Dokl. Akad. Nauk. SSSR, 
320(l):34-39, 1991. 

V.L. Makarov. About functional-discrete approach for solving problems of mathematical 
physics. In Proceedings of Wekua workshop, pages 44-51, Tbilisi, 1993. 

V.L. Makarov. FD-method for nonlinear eigenvalue problems for nonlinear differential 
equations. Dopovidi NAN Ukrainu, 8:16-22, 2008 (ukr.). 

V.L. Makarov and N.. Rossokhata. Error estimates of convergence rate for fd-method for 
sturm-liouville problem with potential in L\. Collected works of Institute of Mathematics 
NAS of Ukraine, 1(3): 1-16, 2005 (ukr). 

V.L. Makarov and N.O. Rossokhata. FD-method for nonlinear eigenvalue problems with 
discontinuous eigenfunctions. Nonlinear Oscillations, 10(1):126-143, 2007. 

V.L. Makarov and N.O. Rossokhata. A rewiev of functional-discrete technique for eigen- 
value problems. Journal of Numerical and Applied Mathematics, 97:97-102, 2009. 

V.L Makarov, N.O. Rossokhata, and B.I. Bandyrskii. Functional-discrete method with a 
high order of accuracy for the eigenvalue transmission problem. Computational Methods 
in Applied Mathematics, 4(3):369-381, 2004. 

V.L Makarov, N.O. Rossokhata, and B.I. Bandyrskii. Functional-discrete method for 
an eigenvalue transmission problem with periodic boundary conditions. Computational 
Methods in Applied Mathematics, 45(2):201-220, 2005. 

V.L. Makarov and O.L. Ukhanev. FD-method for Sturm-Liouville problems. Exponen- 
tial rate of convergence. Applied Mathematics and Informatics (Tbilisi University Press), 
2:1-19, 1997. 



28 



[22] Stephen Nash, David Kahaner, and Clev Moler. Numerical methods and software. 
Prentice-Hall, Inc., New Jersey, 1989. 495 p. 

[23] S.A. Pruess. Estimating the eigenvalues of Sturm- Liouville problems by approximating 
the differential equations. PhD thesis, Purdue University, 1970. 

[24] J.D. Pryce. Numerical solution of Sturm- Liouville problems. Clarendon Press, Oxford, 
New York, Tokyo, 1993. 322P. 

[25] N.O. Rossokhata. Analysis of an eigenvalue transmission problems with FD-method. 
Bulletin of the University of Kiev, (l):194-203, 2006. 

[26] N.O. Rossokhata. FD-method for eigenvalue transmission problems with potential in 
space L\. Bulletin of the University of Kiev, 3:161-168, 2007. 

[27] A. M. Savchuk and A. A. Shkalikov. Sturm-Liouville operators with singular potentials. 
Mat. Zametki, 66(6):897-912, 1999. 

[28] A. M. Savchuk and A. A. Shkalikov. Sturm-Liouville operators with distribution po- 
tentials. Tr. Mosk. Mat. Obs., 64:159-212, 2003. 

[29] Frank Stenger. Handbook of Sine numerical methods. Chapman & Hall/CRC Numerical 
Analysis and Scientific Computing. CRC Press, Boca Raton, FL, 2011. With 1 CD- 
ROM (Windows, Macintosh and UNIX). 

[30] Takayasu MATSUO Tomoaki OKAYAMA and Masaaki SUGIHARA. Error estimates 
with explicit constants for sine approximation, sine quadrature and sine indefinite inte- 
gration. Mathematical Engineering Technical Reports 2009-01, The University of Tokyo, 
January 2009, pages 1-28, 2009. 

[31] Seng V., Abbaoui K., and Cherruault Y. Adomian's polynomials for nonlinear operators. 
Math. Comput. Modelling, 24(l):59-65, 1996. 

[32] V. A. Vinokurov. The eigenvalue and eigenfunction of the Sturm-Liouville problem as 
analytic functions of the integrable potential. Differ. Uravn., 41(6):730-738, 861, 2005. 

[33] V. A. Vinokurov and V. A. Sadovnichii. Asymptotics of arbitrary order of the eigen- 
values and eigenfunctions of the Sturm-Liouville boundary value problem in an interval 
with a summable potential. Izv. Ross. Akad. Nauk Ser. Mat, 64(4):47-108, 2000. 

[34] V. A. Vinokurov and V. A. Sadovnichii. Asymptotics of eigenvalues and eigenfunctions 
and the trace formula for a potential that contains 5-functions. Dokl. Akad. Nauk, 
376(4):445-448, 2001. 

[35] V. A. Vinokurov and V. A. Sadovnichii. Analytic dependence of the eigenvalue and 
eigenfunction of the Sturm-Liouville problem on the integrable potential. Dokl. Akad. 
Nauk, 400(4) :439-443, 2005. 

[36] V.A. Vinokurov and V.A Sadovnichii. The asymptotics of eigenvalues and eigenfunc- 
tions and a trace formula for a potential with delta functions. Differential equations, 
38(6):772-789, 2002. 

[37] V.A. Vinokurov and V.A Sadovnichii. Today theory of the sturm-liouville problem. 
In Proceedings of International Conference 'Tikhonov and contemporary mathematics ', 
pages 383-385, Moscow, 2006. 



29 



